USE  OF  AN  INDEFINITE  INNER  PRODUCT 
FOR  COMPUTING  DAMPED  NATURAL  MODES 


\  by 

B.  N.  Parlett^  and  H.  C.  Chen* 


Abstract  A  quadratic  eigenvalue  problem  with  symmetric  positive  definite  coefficient 
matrices  may  be  reduced  to  linear  form  while  retaining  symmetry  in  the  new  coefficients 
but  neither  of  them  will  be  positive  definite.  Formally  the  symmetric  Lanezos  algorithm 
and  subspace  iteration  may  be  used  to  compute  some  eigenpairs  of  the  linear  problem. 
The  trouble  is  that  the  basis  vectors  are  orthogonal  with  respea  to  an  indefinite  inner  pro¬ 
duct  so  there  is  no  assurance  that  they  arill  be  linearly  independent.  Nevertheless  this  is 
an  attractive  way  to  solve  the  original  problem  and  we  discuss  how  to  implement  it  and 
how  it  relates  to  the  unsymmetric  Lanezos  procedures.  We  discuss  complex  origin  shifts, 
reorihogonalization,  and  error  bounds.  Several  methods  for  solving  the  reduced  problem 
are  mentioned  but  we  have  no  fully  satisfying  technique.  Some  dangers  are  described 
and  examples  are  given  comparing  our  Lanezos  program  with  a  modified  subspace  itera¬ 
tion. 
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1.  INTRODUCTION 


In  the  study  of  Mechanics,  when  one  analyzes  how  small  displacements  from  an  equilibrium  state 
evolve  in  time  one  is  led  to  the  familiar  equations  of  motion 

Mq(r)  +  Cq(/)  +  Kq(r)  =  f(r)  (1.1) 

where  M,  C,  and  K  are,  respectively,  the  /ixn  mass,  damping,  and  stifhiess  matrices  and  q(r),  q(r),  and 
q(r)  are  the  n  xl  acceleration,  velocity,  and  displacement  vectors. 

For  many  solid  structures  attached  to  the  earth  tl^  damping  coefficients  are  small  compared  with 
those  of  mass  and  stiffness.  For  structures  in  space  damping  plays  a  bigger  role  and  the  undamped  model 
may  not  be  a  good  approximation.  To  understand  the  response  of  the  system  to  a  variety  of  external 
forces  f(r)  it  is  still  useful  to  know  the  system’s  dominant  natural  modes  of  vibration.  This  leads  to  the 
quadratic  eigenvalue  problem 

(Xf  M  +  X,  C  +  K)w.  =0  J  =  1 . 2n  (1.2) 

In  general  X,  will  be  complex  and  the  associated  modal  sluq^es  are  given  by  the  real  and  imaginary  parts 
of  w, .  In  the  problems  we  consider  M,  C,  K  are  real,  symmetric,  and  positive  definite. 

Our  paper  has  several  goals.  In  Section  2  we  make  a  plea  for  keeping  symmetry  explicit  in  the 
reduction  to  linear  form  even  though  the  resulting  pencil  (or  pair  of  matrices)  is  not  definite.  We  are  cer¬ 
tainly  not  the  first  to  make  this  suggestion.  In  Section  3  we  show  how  to  make  a  complex  shift  of  origin 
and  yet  keep  the  eigenvalue  algorithm  confined  to  real  arithmetic.  This  idea  appears  not  to  have  been 
used  before.  Sections  4  and  5  are  theoretical,  showing  important  connections  but  do  not  contain  new 
results.  After  that  we  describe  an  implementation  of  a  Lanczos  algorithm  for  symmetric  indefinite  pen¬ 
cils.  It  is  designed  to  exploit  any  small  bandwidth  in  M,  C,  K.  We  describe  how  selective  orthogonaliza- 
tion  can  be  carried  over  from  the  definite  case  (Section  6),  the  numerical  dangers  that  beset  the  algorithm 
in  our  case  (Section  7),  the  solution  of  the  reduced  eigenvalue  problem  (Section  8),  and  how  to  compute 
accurate  error  estimates  at  very  low  cost  at  any  step  of  the  algorithm  (Section  9).  Finally,  some  numerical 
results  and  comparisons  make  up  Section  10. 

We  mention  here  that  certain  problems  need  more  study.  (1)  When  M  is  singular,  or  nearly  so,  the 
computed  Lanczos  vectors  can  acquire  large  components  in  M’s  null  space  without  the  algorithm  being 
able  to  detect  that  this  is  h^pening.  We  have  a  remedy  but  do  not  know  how  completely  it  cures  the 
problem.  (2)  The  efficient  solution  of  the  reduced  eigenvalue  proUem  is  an  interesting  challenge.  How¬ 
ever,  in  the  context  of  the  whole  algorithm,  it  is  satisfactory  to  employ  a  standard,  robust  algorithm  that 
does  not  exploit  aU  the  features  of  the  reduced  problem.  A  Lanczos  program  similar  to  outs  has  been 
developed  independently  (see  [Nour-Omid  &  Regelbrugge,  1988])  and  used  effectively  to  solve  the 
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equations  of  motion  (1.1)  directly  without  computing  eigensolutions. 

We  follow  standard  notational  conventions  in  the  field  of  matrix  computations.  In  particular,  capital 
letters  denote  matrices,  small  Roman  letters  denote  column  vectors,  and  small  Greek  leners  denote 
scalars.  Also  v‘  denotes  the  transpose  of  v  and  a  denotes  the  conjugate  of  a  complex  number  a.  Unless 
the  contrary  is  indicated  llxll  denotes  the  Euclidean  norm  of  x,  ilBII  denotes  the  spectral  norm  of  B. 


2.  REDUCTION  TO  LINEAR  FORM 


There  are  several  ways  of  rewriting  (1.2)  as  a  linear  eigenvalue  problem.  When  K  and  M  are  posi¬ 
tive  definite  then  0  is  not  an  eigenvalue  and  one  reduction  is 


(2.1) 


■  C  M 

'  w,  1  j 

■-K  o' 
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M  0 
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Another  is 


-C  -M 
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K  0 
0  M 


w. 


X,w, 


(2.2) 


The  difference  between  (2.1)  and  (2.2)  may  strike  the  readers  as  frivolous.  But  the  implications  of  this 
change  go  far.  The  reason  that  the  distinction  between  (2.1)  and  (2.2)  seems  negligible  is  that,  at  the  next 
step,  when  these  generalized  problems  are  reduced  to  standard  form  both  (2.1)  and  (2.2)  produce  the  same 
result  since 


(2.3) 
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-1 
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-1 
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This  is  a  nonsymmetric  matrix  as  we  would  expect  from  the  presence  of  complex  eigenvalues. 
Of  course  the  reduction  can  be  made  by  using  the  Cholesky  factorizations 

K  =  Lk  Lr,  M  =  Lm  Lii 


to  produce  first 


Lie*  C  Lit  Lk'  Lm 
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(2.4) 


and  then,  the  standard  form 


-Li?*  C  Li?  -Li?*  L 
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We  suggest  that  this  reduction  to  standaid  fonn,  either  to  (2.3)  or  (2.5),  is  a  tactical  error. 

In  order  to  explain  our  view  the  following  standard  terminology  will  be  needed.  A  symmetric 
matrix  is  indefinite  if  it  has  eigenvalues  of  both  signs. 

Definition  1.  Let  Y  e  be  symmetric  but  indefinite.  The  bilinear  form  defined  by 

( u.  V  )y  :=  v'  Y  u 

is  called  a  pseudo  (or  indefinite  or  improper)  iruier  product  It  obeys  all  the  axioms  of  an  irmer  pro¬ 
duct  except  positivity. 

Definition  2.  Let  F e  R"*^ ,  Y e  R"’* ,  with  Y  symmetric. 

F  is  symmetric  with  respect  to  Y,  if 

F*  Y  =  (YFy=YF  . 


or,  equivalently,  if 


(u,Fv)y  =  (Fu,v)Y 


for  all  u  G  R" ,  V  6  R" , 

In  some  contexts  one  says  that  F  is  self-adjoint  in  Y’s  (pseudo)  inner  product. 
The  trouble  with  (2.2)  and  (2.3)  is  that  neither  of  them  reminds  us  that 


is  symmetric  with  respect  to  both 


-K-‘  C  -K-i  M 
I  0 


C  M 
M  0 


and 


-K  0 
0  M 


Similarly  the  first  matrix  in  (2.5)  is  symmetric  with  respect  to 


but  not  with  respect  to 


-I  0 
0  I 


■  I  0 
0  I 

We  see  no  gain  in  efficiency  in  going  all  the  way  to  (2.5)  in  the  reduction  but  circumstances  may  change. 
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In  brief,  (2.1)  suggests  flie  relevant  symmetries  but  (2.2)  does  not.  (2.1)  is  a  special  case  of  the 
problem 

(A--J-B)z  =  0 

where  A  and  B  are  symmetric  but  indefinite.  A  matrix  pair  (A,  B)  is  sometimes  called  a  matrix  pencil. 
See  [Gantmacher,  1959]  for  example.  In  our  case  we  have  a  symmetric  indefinite  pencil.  Excellent  refer¬ 
ences  for  further  study  of  these  pencils  see  [Gohberg,  Lancaster,  and  Rodman,  1983  and  1986]. 

Lemma.  B“*  A  is  symmetric  with  respect  to  both  A  and  B.  A  B“*  is  symmetric  with  respect  to 

both  A"’  and  B~’. 

Proof:  For  all  u  ,  v  in  IR" , 

(  u,  B~*  A  V  )a  =  v*  A  B"*  A  u  =  (  B"’  A  u,  v  )a  . 

(  u,  B“*  A  V  )b  =  v*  A  u  =  (  B“*  A  u,  V  )b  . 

The  second  result  can  be  obtained  similarly.  □ 

It  turns  out  that  the  Lanczos  algorithm  may  be  invoked  with  the  operator  B~*  A  using  a  pseudo  inner 
product  defined  by  A  or  B.  The  three  term  recurrence  still  holds  in  this  more  general  situation.  There  is 
more  on  the  Lanczos  algorithm  in  Section  4. 

Subspace  Iteration  [Bathe  and  Wilson,  1976]  or  [Rutishauser,  1970]  (also  called  the  method  of 
Simultaneous  Iterations)  also  extends  formally  to  the  indefinite  case  but  the  Rayleigh  Ritz  approximations 
produced  at  each  step  are  no  longer  optimal  in  any  meaningful  sense.  We  have  made  the  necessary 
modification  to  our  standard  SI  program  and  use  it  as  a  simple  rival  to  our  Lanczos  procedure. 

In  contrast  to  the  definite  case  both  algorithms  can  break  down  or  become  unstable  when  close  to 
breakdown.  Alarming  things  can  happen  with  an  indefinite  (or  improper)  inner  produa  :  a  set  of  orthogo¬ 
nal  vectors  might  be  linearly  dependent.  The  geometry  associated  with  such  an  improper  inner  product  is 
the  geometry  of  relativity  theory. 

Next  we  wish  to  point  out  that  it  is  possible  to  work  with  A  B"’  instead  of  B~’  A  but  we  see  no 
advantage  to  this  formulation  in  our  problem.  Since  A  B~'  is  symmetric  with  respect  to  B~'  and  A"'  we 
must  work  with  these  pseudo  inner  prxxlucts.  Now  A"*  is  complicated  but  when  B  =  diag  (-K,  M)  then 
we  can  form  the  Cholesky  factorization  of  M  and  K  once  for  all.  After  that  the  computation  of  B"'  v 
costs  no  more  than  the  computation  of  B  v.  The  need  for  shifts  of  origin  (see  next  section)  suggests  the 
use  of  an  operator  A  and  the  pseudo  irmer  product  defined  by  A.  Here  H  may  be  B  or  B  -  o  A  or 
Re  (B  -  o  A'/  or  Im  (B  -  o  A). 
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Both  fonns  (2.1)  and  (2.2)  are  well  known.  In  a  recent  paper  [Borri  &  Mantegazza  1977]  fonn  (2.1) 
is  explicitly  proposed  when  M,  C,  and  K  are  all  symmetric.  However,  no  mention  is  made  of  the  fact  that 
the  Rayleigh  quotients  used  in  that  paper  may  overflow  or  yield  0/0.  A  conventional  alternative  to  the 
pseudo-symmetric  form  (2.1)  or  (2.4)  is  to  apply  the  two  sided  Lanczos  algorithm  to  the  matrix  in  (2.3)  or 
(2.5).  In  Section  4  we  show  the  connection  between  these  two  approaches. 

We  recall  that  the  attraction  of  the  Lanczos  algorithm  is  that  a  few  of  the  largest  eigenvalues  may  be 
found  by  stopping  the  algorithm  long  before  the  lull  n  steps.  It  is  usual  practice  to  give  the  algorithm  an 
operator  such  as  A  whose  largest  eigenvalues  are  shifted  reciprocals  of  the  ones  we  really  want.  We 
can  expect  to  find  the  p  largest  eigenvalues  within  max  {p+8, 2p }  steps.  See  Section  10  for  examples. 
By  reducing  the  n  xn  quadratic  problem  to  linear  fonn  the  dimension  of  the  Lanczos  vectors  becomes  2n . 
However  the  cost  of  implementing  the  Lanczos  algorithm  is  only  doubled  (approximately)  because  the 
structure  of  the  matrices  A  and  B  may  be  exploited.  The  algorithm  and  operation  counts  are  given  in  Fig¬ 
ure  1. 

Note  that  we  could  have  used  the  algorithm  for  the  operator  H~’  A  using  the  pseudo-inner  product 
defined  by  H.  The  variation  alters  the  subroutines  supplied  to  the  Lanczos  algorithm,  not  the  program 
itself.  There  may  be  advantages  to  the  H  inner  produa  but  we  have  not  studied  the  matter  since  H  may 
be  complicated. 


3.  ORIGIN  SHIFTS 

There  are  applications  (e.g.  space  structures)  in  which  the  stiffiiess  matrix  K  is  singular.  This 
makes  B  =  diag  (-K,  M)  singular  too.  In  this  situation,  we  need  to  solve  a  shifted  problem 

(X-o)  Az  =  (B-oA)z 

where  o  ^  0  is  a  real  shift.  To  preserve  the  block  diagonal  form,  we  note  that  the  B  -  o  A  can  be  factored 
into 

■  I  -ol 
0  I 

We  need  to  factor  only  M  and  the  shifted  stiffness  matrix  K  o  C  M,  which  has  the  same  banded 
structure  as  K  in  order  to  solve  (B  -  a  A)u  =  Av  for  u  given  v.  Thus  we  may  solve 

(B  -  oA)"'Az  =  — -z  to  obtain  eigenvalues  close  to  o.  This  is  standard  practice. 

(X-o) 


(-K-oC-o^M)  O] 

I  0 

0  Ml 

-ol  I 

In  several  applications  the  user  would  like  to  explore  pan  of  the  complex  plane  and  see  how  cenain 
eigenvalues  vary  with  changes  in  the  matrix  elements.  Unfortunately,  if  o  is  complex  then  B  -  o  A  is  not 
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real  and  so  (B  -  o  A)"’  A  is  not  a  valid  operator  for  our  Lanczos  algorithm. 

A  way  out  of  this  quandary  was  proposed  in  [Pailett  and  Saad.  1987].  The  operator  given  to  the 
Lanczos  algorithm  is  either 

Re[(B-oAr']A  or  Im  [  (B  -  o  A)"’ ]  A  .  (3.2) 

Here  Re  G  (G  +  G),  Im  G  :=  —7;-  (G  -  G),  and  G  is  the  conjugate  of  G.  Both  operators  are  real 
2  2i 

and  have  complex  conjugate  pairs  of  eigenvalues. 

Let  z  be  an  eigenpair  of  (B,  A),  i.e.  they  satisfy  A.  A  z  =  B  z.  By  subtracting  a  A  z  and  a  A  z 
from  each  side  and  inverting  one  finds  that 

Re  [  (B  -o  A)-'  1  A  z=  I  +  7^)  z  .  (3.3a) 

2  X-a  X-a 

Im[(B-oA)-MAz  =  ;^(Y^-7^)z  .  (3.3b) 

2 1  A-<j  X-a 

To  recover  X  from  the  computed  eigenvalues  of  either  of  these  operators  requires  merely  solving  a  qua¬ 
dratic  equation  for  the  root  that  is  closer  to  a.  However  X  may  be  recovered  directly,  with  extra  work., 
from  the  Rayleigh  quotient 

Xs?  Bz/z*  Az  . 

In  order  to  exploit  narrow  bandwidth  it  is  convenient  to  use  complex  arithmetic  in  (3.1)  to  factor 
B  -  o  A  and  thus  to  solve  ,  for  w  e  C  " ,  (B  -  o  A)  w  =  A  v  for  any  given  v  e  IR" .  However,  the  subrou¬ 
tine  returns  either  Rew  or  Imw  and  so  there  is  no  complex  arithmetic  in  the  Lanczos  algorithm.  The 
major  computation  is  in  factoring  K +0  C  M.  This  is  done  once.  Note  that  the  pseudo-inner  product 
is  defined  by  A  and  is  independent  of  a. 

Both  the  real  and  imaginary  part  of  (B  -  o  A)"‘  provide  good  enhancement  of  the  eigenvalues  close 
to  a.  It  is  shown  in  [Parlett  and  Saad,  1987]  that  the  imaginary  part  suppresses  the  unwanted  eigenvalues 
far  from  o  more  strongly  than  does  the  real  part  More  experience  is  needed  on  this  aspect  of  the  shift. 

4.  CONNECTION  BETWEEN  LANCZOS  ALGORITHMS 

The  original  formulations  given  by  Lanczos  (in  [Lanczos,  1950])  considered  the  standard  symmetric 
eigenvalue  problem  A  z  =  X z,  i.  e.  B  =  I.  The  Lanczos  vectors  {q*  )t,i  for  symmetric  nxn  A  obey  the 
well-known  3-teTm  recurrence  relation 


A  *1*  =  PtQjk-i  +  tt*q*  +  Pt>iqii4-i 


/fc  =  1, . . . ,  n 


(4.1) 
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with  qo^O.  The  {qtijt.i  form  an  onhononnal  set.  The  attractive  feature  of  the  Lanczos  algorithm  for 
the  generalized  problem  A  -  X  B,  where  B  is  symmetric,  positive  definite  is  that,  as  in  (4.1) 

A  q*  =  P*<U-i  +  +  Pt+iQi+i 

The  Lanczos  vectors  are  still  oithonormal  but  in  the  inner  product  defined  by  B.  Usually  they  will  not  be 
orthogonal  in  the  ordinary  Euclidean  sense. 

In  general,  for  a  nonsytmnetric  matrix  F  (we  rename  B~'  A  as  F)  there  is  a  two-sided  generalization 
(also  given  by  Lanczos)  in  which  two  sets  ( Pt, .  - . .  Pm  )  and  { qj, ....  qM  }  are  generated.  Moreover,  if 
pi  <Lk  =  f^cn 

0*  =  P**!*-!  +  +  Yk+ltU+l  • 

F'  P*  =  Y*P*-i  +  a*P*  +  P*+iP*+i .  (4.2) 

and  the  { P;t  } ,  ( q^  }  are  biorthogonal, 

p/q*=0.  i*k 

This  algorithm  will  break  dovm  when  pi  q^  =  0  with  pj^  ^  0,  q^  #  0.  We  can  expect  inaccuracy  when 
p*  q*  is  very  small.  There  is  more  on  this  point  in  Section  5. 

For  theoretical  purposes  it  is  convenient  to  normalize  the  Lanczos  vectors  by 

Piq*  =  l  llp*ll  =  llq;kll  (4.3) 

provided  breakdown  does  not  occur.  However,  in  practice,  it  is  better  to  keep 

||p^ll  =  11(^11  =  1 

and  define 


®*=piq* 


(4.4) 


If  to*  is  tiny  (  <  KT*  or  e*'^)  then  the  { p*  }  and  { q*  }  are  nearly  linearly  dependent  and  consequently 
make  bad  bases  in  which  to  represent  the  solutions.  The  normalization  (4.3)  conceals  any  deterioration  in 

the  quality  of  the  basis  { qj . qM  ).  The  normalization  (4.4)  makes  the  three  term  recurrence  slightly 

more  complicated.  Define 


Yz 

p2  02  Yj 


Jm  = 


Pm-J  O,.-!  Ym 
Pm  Om 


=  PmFQm. 
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Q, -diagimi . Q».. 

Q-  =(qi.  -  .,  <1.). 

P«=(Pi . P«)  (4.5) 

then,  in  matrix  form,  the  three  term  recurrence  is 

FQ,  -Q,  Q;*  J,  =0  , 

Pir>j,  Q;»Pi=o.  (4.6) 

Note  that  when  F  =  A  =  A'  and  B  =  I  then  P„  =  Q, ,  =  I,  and  J„  =  T„  =  . 

In  our  pseudo-symmetric  eigenvalue  problem  A  -  XB  we  face  the  same  possibility  as  in  (4.2)  that 
the  Lanczos  process  may  break  down  and  may  be  unstable  when  close  to  breakdown.  To  distinguish  this 
algorithm  from  the  two-sided  algorithm  given  by  (4.5)  and  (4.6)  we  label  the  Lanczos  vectors  in  our  algo¬ 
rithm  as  qi, . . . ,  ^ .  It  is  important  to  know  when  this  instability  occurs  and  so  we  prefer  to  normalize 
the  qjk  to  satisfy 

llq*ll  =  l  ,  for  aU  k  ,  (4.7) 

and  then  define 

=(4k»  4k)A  ■ 

A  tiny  value  of  ©*  is  a  sign  of  danger.  The  three  term  recurrence,  in  matrix  form,  is 

AQ,  -Q,  =0  , 

f,  =  A  B-*  A  Q„  ,  (4.8) 

Q„=Q‘AQn 

provided  that  breakdown  does  not  occur.  Note  that  t„  is  symmetric.  Equating  the  kth  column  on  each 
side  of  (4.8)  3delds  the  three  term  recurrence  for  our  pseudo-symmetric  formulation 

B~‘  A  q*  =  (0*/(n*_i)  4k-i  +  (o*/<Ojk)  q*  +  (0*+i/to*+i)  ,  k  <  n  .  (4.9) 

This  equation  shows  the  possible  danger  of  small  values  among  the  {m, }. 

The  goal  of  the  preceding  discussion  was  to  prepare  for  the  interesting,  but  natural,  result  that  our 
pseudo-symmetric  procedure  (4.8)  &  (4.9)  is  a  disguised  form  of  the  two-sided  algorithm  (4.6).  By 
choosing  Pi  appropriately  it  turns  out  that  p^  =  A  q^^,  for  all  k,  and  so  there  is  no  need  to  hold  the  p^ 
explicitly  in  memory. 
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Theorem.  Suppose  that  A  and  B  are  real  and  symmetric.  Suppose  that  no  breakdown  occurs.  Let 
the  two-sided  Lanczos  algorithm  in  (4.6)  be  applied  with  operator  B~'  A  (starting  with  pi  and  q^)  to  pro¬ 
duce  two  bioithogonal  sequences  {ptlt.i,  {qAlt.i.  Let  the  pseudo-symmetric  Lanczos  procedure  in 
(4.9)  be  applied  with  B”*  A  (starting  with  qj)  with  pseudo-inner  product  defined  by  A  to  produce  {q*  )*_!. 
If  qi  qi  and  pj  =  A  qj  /  IIA  qill  then  4k  =  q^  and  pj^  -  A  /  IIA  qj^ll,  it  =  2 . n. 

Proof  :  We  use  induction  and  let  u  I  I  v  mean  that  u  is  a  nonzero  multiple  of  v.  Put  it  =:  1  in  (4.9)  and 
equate  column  1  on  each  side  of  (4.6)  to  find 

q2  I  I  B"‘ Aqi--qi(di/co,)  , 

q2  I  I  B~' Aqi-q,  (tti/coi)  , 

P2  I  I  A  B“’  p,  -  pj  (ai/to,) . 

By  choice  q^  =  qj  and  A  qj  =  pj  IIA  qill,  so 

:=  (Qi.  AOa  =  (qt.  qi)A  =  q{  Pi  "A  qiH  =  ©1  IIA  q,ll  . 

dj  ;=  (qi,  B"‘  A  q^^  =  (qj.  B"'  A  qi)^  =  Pi  B"'  A  qj  IIA  qjll  =  a,  IIA  q,ll . 

Thus  di/(Oi  =  tti/Wj  and  hence  q2  =  qj  and  pj  1 1  A  q2  /  IIA  qjll.  So  pj  =  A  q2  /  IIA  q2ll. 


The  induction  assumption  is  that  the  theorem  holds  for  /fc  =y  -  I  and  j.  Put  k  -  j  in  (4.9)  and 
equate  column  J  on  each  side  of  (4.6)  to  find 


4+1  I  I  B-'  Aqy  -qj  (d/(0;)-qy.,  (^/w^.i)  . 
qy+,  I  I  B-'  Aq^  -qj  (a/wp-qy.,  (y,/©^.,)  . 
Py+i  I  I  AB"'  p^  -pj  (ay/©;)-p^_i  (P>/©y_,) . 
Use  the  induction  assumption  to  verify  that 


(0;_i  :=  (qy-i.  4-i)a  =  (qy-i.  q;-i)A  =  qj-i  P;-i  HA  q;_ill  =  IIA  qy_ill  , 

%  :=  (q;-i-  ^  q>)A  =  (q>-i.  B-‘  A  q^A  =  P/-i  B"'  A  q^  IIA  qy_,ll  =  IIA  q^.j! 

:=  (qy  -  4 Oa  =  (qy •  q>)A  =  qJ  Py  ha  qyll  =  ©;  HA  q^ II  . 
tty  :=  (qy .  B"’  A  q^A  =  (qy .  B"‘  A  q;)A  =  p]  B’’  A  q^  IIA  q^ II  =  IIA  q^  II  , 

_1IA^  , 


,,  :=  p;.,  B-'  A  „  =  IIA  q,.,ll  A  B-  A  ,,  =  q'.,  A  B-  p,  =  P, 


using  these  relations  we  verify  that 


qy>J”qy+i  ’ 
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Py+1  I  ^  ^  fly+1  • 

Thus,  =  A  qj^i  /  IIA  qj^iW  and  the  result  holds  for  k-j+\.  By  the  principle  of  induction,  if  the 
algorithms  do  not  break  down  then  the  theorem  holds  for  k  =  1. 2 . n.  □ 

5.  MOMENT  MATRICES  AND  HYPERBOUC  PAIRS 

The  breakdown  of  the  (generalized)  form  of  the  Lanczos  algorithm  may  be  interpreted  as  break¬ 
down  in  the  triangular  factorization  of  a  certain  matrix.  This  connection  is  well  known  (see  [House¬ 
holder,  1964]  or  [Parlett,  Taylor,  and  Liu,  1985])  and  gives  valuable  insight  We  review  it  now.  Let 

K,  =[q„B-‘Aqi,(B-'A)2qi . (B-^Ay-%]. 

This  is  called  a  Krylov  matrix.  When  successful  the  Lanczos  algorithm  constructs  a  matrix 

Q>  =  [qi.  -.  q>] 

of  Lanczos  vectors  which  are  A-orthogonal,  i.e. 

Qj  A  Qj  =  Cl j=  diagonal. 

Moreover  Qy  and  Ky  are  connected  by 

Ky=QyL; 

where  Ly  is  an  invertible  lower  triangular  matrix.  In  other  words  Qy  is  obtained  from  Ky  by  the  Gram- 
Schmidt  process.  However,  Lanczos  found  a  clever  way  to  find  qy  without  actually  invoking  Gram- 
Schmidt.  This  result  is  not  obvious:  it  expresses  the  fact  that  B"'  A  q^,,  is  a  linear  combination  of  q^_i, 
q„, ,  and  q^+j  and  from  this  it  can  be  deduced  that  for  each  m  =  1, 2 . j,  (B~'  A)”  qj  is  a  linear  com¬ 
bination  of  qj,  q2 . q„,+i  only. 

In  the  present  context  the  moment  matrix  is  defined  as 

My  =[mj*]  , 

mu,  =  (q{,  (B-‘A)'^*-2  q,)A  =  q{  A  (B-'A)'^*-^  q,. 

A  little  manipulation  shows  that 

My=K;AKy  . 

Consequently,  if  the  Lanczos  process  does  not  breakdown. 


M^  =  Ly  q;  a  Qy  l; = hj  Qj  l;  . 

This  is  the  triangular  factorization  of  M^.  We  do  not  insist  that  the  diagonal  elements  of  L^  be  1. 
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Triangular  factorization  exists  if  all  the  leading  principal  submatrices  of  are  nonsingular  (i.e.  inverti¬ 
ble).  Cbnversely,  if  the  last  diagonal  element  of  Qj  is  the  first  to  vanish  then  the  Lanczos  algorithm 
breaks  down  at  the  end  of  step  j . 

Observe  that  Kj  and  are  each  nx;  matrices  whereas  Lj  is  invertible.  Thus  has  full  rank  if 
and  only  if  Qj  has  full  rank.  It  is  desirable  (though  not  essential)  that  Qj  have  full  rank. 

Lemma.  Let  ,  Mj ,  and  A  be  as  defined  above.  If  Ky  has  lull  rank  and  A  is  symmetric,  posi¬ 
tive  definite  then  so  is  and  the  factorization  Lj  Qj  Lj  exists  with  Qj  positive  definite. 

Proof :  For  any  v  e  IR-' . 

v*  My  V  =  v'  Ky  A  V  =  w*  A  w  >  0  , 

unless  w  =  0.  If  Ky  has  full  rank  then  w  s  0  implies  v  =  0.  Hence  My  is  positive  definite.  The  leading 
principal  submatrices  of  a  positive  definite  matrix  are  positive  definite  and  thus  My  permits  a  triangular 
factorization.  □ 

When  A  is  indefinite  then  there  exist  veaors  w  «  0  such  that  w'  A  w  ^  0  and,  except  in  trivial  cases, 
there  win  be  starting  vectors  qj  and  j -vectors  v  such  that  the  A-norm  of  Ky  v  vanishes  for  large  enough 

j- 

Definition:  v  is  isotropic  if  v'  A  v  =  (v,  v)^  =  0. 

It  is  helpful  to  regard  the  occurrence  of  breakdown  (i.e.  when  qy  is  isotropic)  not  as  disaster  but  as  a 
reminder  that  A  is  indefinite.  There  is  a  natural  way  to  proceed.  The  idea  goes  back  to  Lagrange  but  was 
pui  to  use  by  D.  G.  Luenberger  to  extend  the  conjugate  gradient  algorithm  [Luenberger,  1969]  and  by 
[Bunch,  Kauffman,  and  Pariett,  1976]  to  stablize  triangular  factorizatioa 

By  a  simple  rotation  of  coordinates  the  hyperbola  -  y  ^  =  1  may  be  written  as  ^  q  =  2.  In  our  con¬ 
text  the  idea  is  to  modify  our  generalized  Lanczos  algorithm  as  follows.  If  qy  is  an  isotropic  vector  then 
choose  the  value  of  in  (4.1)  or  (4.9)  so  that  qy+j  is  also  isotropic.  Provided  that  qy  A  qy+i  is  not  too 
small  these  Lanczos  vectors  still  provide  a  good  basis.  They  are  no  longer  A-orthogonal.  Nevertheless 
Qy'+i  A  Qy+i  is  block  diagonal  with  2x2  diagonal  blocks  of  the  form 

'  0  ©yy,,- 

0 

k  s 

whenever  the  modification  is  used.  We  follow  Luenberger  [Luenberger,  1969]  in  calling  qy,  qy+i  a 
hyperbolic  pair. 
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In  practice  there  is  no  need  to  wait  until  qj  A  is  negligible  before  switching  to  a  hyperbolic  pair 
for  qy  and  qy>i.  By  forsaking  strict  A-orthogonality  we  gain  a  better  ctmditioned  basis.  Similar  thinking 
for  the  two  sided  algorithm  yields  die  Look  Ahead  algorldim  [Parlett,  Taylor,  and  Liu,  1985]  but  the  pro¬ 
cedure  is  ladier  complicated.  Our  experience  is  limited  but  the  only  breakdowns  we  have  encountered 
were  produced  deliberately. 


6.  LOSSOFORTHOGONAUTY 

It  is  well  known  that  the  three  term  recurrence  does  not  produce  A-oithogonal  Lanczos  veaors  in 
finite  precision  arithmetic.  It  is  the  convergence  of  the  algorithm  that  provokes  this  deterioration  not  can¬ 
cellation.  There  are  three  ways  to  respond  to  the  situation. 

1.  Do  nothing. 

This  technique  causes  Lanczos  sequences  to  be  2  or  3  times  larger  than  necessary.  It  is  only  of 
interest  when  the  whole  spectrum  is  wanted.  The  loss  of  orthogonality  does  not  prevent  the  calculation  of 
fully  accurate  eigenvalues  and  eigenvectors.  It  merely  slows  down  the  process. 

2.  Full  reorthogonalization  at  each  step. 

Here  the  vector  computed  by  the  three  term  recurrence  is  explicidy  orthogonalized  against  all 
preceding  Lanczos  vectors.  This  requires  keeping  the  auxiliary  vectors  (=  Aq*)  in  fast  memory 
(unless  virtual  memory  is  in  use)  and  using  them  every  step.  However  for  short  runs  of  20  or  30  Lanczos 
steps  the  cost  is  not  excessive. 

3.  Selective  reorthogonalization. 

Most  of  the  benefits  derived  from  a  set  of  Lanczos  vectors  that  is  A-orthogonal  to  working  accuracy 
are  also  enjoyed  when  the  vectors  are  merely  semi-orthogonal. 

Definition.  Let  e  be  the  roundoff  unit.  Then  two  vectors  u,  v  are  semi-orthogonal  with  respect  to  A 
if 

1  u*  A  v  I  <  Ve  Hull  •  iivii  •  HAM  . 

Two  techniques  have  been  proposed  for  maintaining  semi-orthogonality.  In  [Paiiett  and  Scott,  1979] 
it  is  shown  that  computed  vectors  tend  to  be  pulled  towards  converged  Ritz  vectors.  HetKe  it  helps  to 
orthogonalize  the  current  Lanczos  vectors  against  such  vectors  from  time  to  time.  This  remedy  only 
addresses  one  mechanism  that  provokes  orthogonality  loss. 
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In  [Simon,  1984]  a  recunence  was  found  that  governs  the  inner  products  among  the  Lanczos  vectors 
themselves.  This  recurrence  permits  an  accurate  estimate  of  the  orthogonality  loss  in  the  current  Lanczos 
vector  to  be  computed  at  each  step  at  a  cost  proportional  to  the  number  of  Lanczos  steps  taken.  It  is  then 
easy  to  orthogonalize  the  new  vector  against  all  the  old  ones  when,  and  only  when,  necessary  in  order  to 
maintain  semi-orthogonality. 

The  recurrence  extends  easily  to  cover  A-otthogonality  and  our  algorithm  incorporates  this  form  of 
selective  orthogonalizatioa  As  a  rule  of  thumb  selective  orthogonalization  has  a  cost  about  1/3  the  cost 
of  full  reorthogonalizadon.  For  large  problems  on  some  computer  systems  the  I/O  costs  of  reorthogonali- 
zation  dominate  the  arithmetic  costs. 

A  full  discussion  of  semi-orthogonality  and  how  to  maintain  it  is  given  in  [Pailett,  Nour-Omid,  and 
Liu.  1988], 


7.  DANGERS 

The  previous  sections  show  that  the  symmetric  quadratic  eigenvalue  problem 

(A,^M4-A,C  +  K)u  =  0  (5.1) 

may  be  reduced  to  the  simpler  form 

(T-^Q)u  =  0  (5.2) 

where  T  is  real  symmetric  and  tridiagonal  and  Q  is  real  diagonal  (or  block  diagonal). 

Any  difficulties  inherent  in  the  original  problem  must  be  inherited  by  the  reduced  problem.  Impor¬ 
tant  feamres  are : 

(1)  complex  eigenvalues  are  present  and  are  wanted. 

(2)  Sometimes  (depending  on  C)  an  eigenvalues  X  may  be  degenerate;  it  may  belong  to  a  nondiagonal 

Jordan  block  of  T.  This  can  happen  when  X.  is  a  double  real  eigenvalue  (the  coalescence  of  a 
conjugate  pair  of  eigenvalues)  or  when  both  X  and  X  are  double  eigenvalues. 

We  consider  (1).  Since  T  and  Q  are  real  it  is  desirable  to  postpone  the  use  of  complex  arithmetic. 
Indeed  if  the  HR  algorithm  [Bunse-Gerstner,  1981;  Brebner  and  Grad,  1982]  is  used  then  the  pair  (T,  Q) 
is  eventually  transformed  into  ( t,  A )  where 

A  =  diagi±\,±\ . ±1) 

and  t  is  block  diagonal  with  2x2  and  1x1  blocks.  Each  complex  conjugate  pairs  of  eigenvalues  is  found 
from  a  real  pencil  of  the  form 


( 


No  complex  arithmetic  is  needed  in  this  case  but  the  HR  algorithm  is  not  always  stable.  However  there 
are  alternative  techniques  for  finding  some  or  all  of  the  eigenvalues  of  (T,  Q)  and  exploiting  the  banded 
form.  SeeSecUonS. 

Thus  (1)  is  not  a  serious  difficulty. 

Now  consider  (2) :  degenerate  multiple  eigenvalues.  Theorem  lS-2-1  in  [Parlett,  1980]  states  that 
any  real  square  matrix  may  be  written  in  the  form  A  where  A  and  B  are  real  symmetric.  Conse¬ 
quently  a  symmetric  indefinite  pair  (A,  B)  may  suffer  from  highly  defective  eigenvalues.  A  simple  exam¬ 
ple  is 


'  0  0  1  a' 

'  0  0  0  r 

0  1  a  0 

0  0  10 

1  a  0  0 

-X 

0  10  0 

a  0  0  0 

.  j 

10  0  0 

•.  1 

which  is  a  Jordan  block  in  disguise. 

If  T  has  some  off  diagonal  elements  that  vanish  then  the  reduced  eigenvalue  problem  splits  up  into 
smaller  pieces.  Each  piece  consists  of  a  T  with  all  its  values  nonzero. 

Lemma:  If  no  off  diagonal  entry  P,-  of  T  vanishes  then  to  each  eigenvalue  of  T  -  (X  -  o)“*  Q  there 
corresponds  exactly  one  eigendirectioa 

Proof:  The  minor  of  the  (i ,  j)  entry  is  (pjPj  -Py-i)  ^  0  for  all  X  and  so  the  rank  never  drops  below  y  -1 . 

□ 

However  it  is  still  possible  to  have  eigenvalues  of  high  multiplicity. 

This  result  is  to  be  expected.  Whatever  the  geometric  multiplicity  of  an  eigenvalue  of  A  -  XB  the 
Lanczos  algorithm  can  only  "see"  the  projection  of  the  starting  vector  qj  onto  the  invariant  subspace  asso¬ 
ciated  with  the  eigenvalue.  Thus,  in  exact  arithmetic,  it  is  possible  for  T  -  X  Q  to  have  generalized  eigen¬ 
vectors  with  lower  grades  than  the  tme  multiplicity  of  the  associated  eigenvalues.  It  is  exactly  the  same 
in  the  truly  symmetric  case.  The  Lanczos  algorithm  cannot  "see"  geometric  multiplicities. 

Unfortunately  it  is  the  prescence  of  degenerate  eigenvalues  that  causes  breakdown  of  the  HR  algo¬ 
rithm. 


Example  of  a  multiple  eigenvalue  : 
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has  0  as  a  double  eigenvalue  with  eigenvector  (1.-1)'  and  any  vector  other  than  the  eigenveaor  is  a  gen¬ 
eralized  eigenvector  of  grade  2.  There  is  no  simpler  symmetric  representation  of  such  an  eigenvalue. 

We  are  studying  the  effect  of  multiple  eigenvalues  on  the  performance  of  the  Lanczos  algorithm. 

There  is  another  difficulty  more  insidious  than  the  first  two. 

(3)  Undetected  growth  of  Lanczos  vectors  in  certain  directions.  The  following  example  was  given  by  Dr. 
T.  Ericsson  [private  coimnunicadon] 

Example : 

A.  =  diag{\,-\.x,x . x), 

B=diag(-l,  l.x.x, ,. . . ,  x), 
qi  =  (o,o.x.x . x)‘. 

Then  (q^  qi)^  is  independent  of  o  and  cotvsequently  both  T  and  Q  will  be  independent  of  a.  The  trouble 
is  that 

a =(1,1,0 . oy 

is  an  eigenveaor  of  B~*  A  with  eigenvalue  I  and  (z,z)A  =  (a.  a)B  =  0.  Arbitrarily  large  multiples  of  z 
could  be  present  in  the  Lanczos  vectors  (on  account  of  qj)  and  the  Lanczos  algorithm  would  be  blind  to 
them.  This  possibility  of  undetectable  growth  in  certain  directions  is  a  generalization  of  the  phenomenon 
reported  in  [  Nour-Omid,  Paiiett,  Ericsson,  &  Jensen,  1987]  where  the  direction  was  in  the  nuD  space  of 
A.  Here  it  is  the  isotropic  eigenvectors  that  are  invisible. 

It  remains  to  be  seen  whether  the  practice  of  keeping  all  Lanczos  veaors  with  Euclidean  length  1 
will  alleviate  the  problem  or  merely  drive  the  Lanczos  vectors  into  the  space  panned  by  undetectable 
eigenveaors. 


8.  HOW  TO  SOLVE  T  -  \  Q 

(1)  Use  EISPACK  on  O"*  Tj.  (Subroutines  HQRl  or HQR2) 

The  tridiagonal  form  will  be  expanded  to  Hessenberg  form  by  the  QR  algorithm.  Thus  no  advan¬ 
tage  is  taken  of  the  compact  tridiagonal  form.  The  arithmetic  effort  is  approximately  to  find  all  the 
eigenvalues  at  step  j. 

(2)  The  HR  algorithm  [Bunse-Gerstner,  1981;  Brebner  and  Grad,  1982). 
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Let  =  1(0, 1,  1=1,...,;,  and  A  =  diVz^(Si . 5j).  The  pencil  (T, Q)  is  equivalent  to 

(A~*  T  A"*,  i)  where  1^,  =  eOj  /  I  eo,  I ,  i  =  1, . . . ,  ; .  TTie  HR  algorithm  is  a  generalization  of  the  sym¬ 
metric  tridiagonal  QR  algorithm.  However  hyperbolic  rotations  of  the  form 


cosh0  sinh0 

sec  0  tan  0 

sinh0  cosh0 

or 

tan  0  sec  0 

are  used  in  place  of  plane  rotations  wheiMver  the  associated  diagonal  elements  of  I  have  opposite  signs. 
Complex  shifts  may  be  used  without  evoking  complex  arithmetic  in  the  same  way  as  they  are  used  in  the 
unsymmetric  QR  algorithm. 

The  only  weakness  is  that  the  HR  algorithm  can  break  down  for  certain  shift  values  and  can  also  be 
unstable  when  the  shift  is  close  to  breakdown.  In  addition  it  will  find  all  the  eigenvalues  at  each  step. 
This  is  overkill.  When  the  algorithm  succeeds  it  needs  aj^roximately  10;^  arithmetic  operations  for  a 
JxJ  pair  (T.  ft). 

(3)  The  LR  algorithm  applied  to  Q"*  Ty. 

This  is  closely  related  to  the  HR  algorithm.  It  is  efficient,  finds  all  the  eigenvalues  and  can  break 
down. 

(4)  analyze; 

This  application  of  Laguerre’s  algorithm  was  developed  for  use  with  the  two-sided  Lanczos  algo¬ 
rithm.  A  data  structure  consisting  of  some  eigenvalues  of  the  pencil  (Ty_p  Qy-i)  together  with  their  error 
estimates  is  updated  to  produce  some  eigenvalues  of  (Ty,  Qy)  and  their  error  estimates.  These  estimates 
are  measures  of  how  close  the  eigenvalue  is  to  an  eigenvalue  of  the  origiiral  pair  (A,  B-oA). 

To  update  one  eigenvalue  a  sequence  of  Laguerre  iterates  is  computed  starting  with  an  old  one.  The 
calculation  of  the  Laguerre  iterate  takes  lull  advantage  of  the  tridiagonal  fonn  but  complex  arithmetic  is 
used.  The  program  is  still  under  development.  The  difficulty  is  to  detea  and  find  new  large  eigenvalues 
that  are  not  close  to  any  at  the  previous  Lanczos  step.  One  remedy  is  simply  to  find  all  eigenvalues  of 
(Ty ,  Qy )  at  each  step. 


Let 


9  ERROR  BOUNDS 

(Ty  -0Qy)S  =  O  . 


(9.1) 


with  llsll  =  1 ,  define  a  typical  eigenpair  (0,  s)  of  the  reduced  problem.  Let  F  be  the  operator  given  to  the 
Lanczos  algorithm  and  let  F  be  symmetric  with  respea  to  A.  After  ;  steps,  in  exact  arithmetic. 
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F  Qj  -  Qj  Q,-‘  T,  =  (P^^,  /  <Oj^0  ej  (9.2) 

The  ‘Ritz  vectors’  corre^nding  to  6  in  (9.1)  is  defined  by 

y:=QyS.  (9.3) 

To  obtain  an  enor  bound  for  the  approximate  eigenpair  (6,  y)  one  postmultiplies  (9.2)  by  s  and  uses  (9.1) 
to  obtain 

F  y  -  y  0  =  qy+,  (py^,  /  (Oy^,)  s0‘ )  (9.4) 

and 

IIF  y  -  y  011  =  (py.,1  /  <0y+,)  I  sO)  I  (9.5) 

since  llqy+ill=  1.  Unfortunately,  it  is  not  possible  to  evaluate  llyll  without  computing  y.  Note  that 
f  A  y = ^  Qj  s  but  this  quantity  might  be  0  or  negative  in  general. 

When  A  is  positive  definite  then  it  defines  a  proper  norm  on  R"  via 

IIvIIa  =  (v'  a  v)*'^ 

InthatcasellqillA  =  V^,  j  =1, ....  j,  and 

IIF  y  -  y  0iIa  =  (P;>i '  I  *0)  I  . 

llyllA  =  (s'Qys)‘^  .  (9.6) 

The  following  error  bound  extends  well  known  results  for  symmetric  matrices.  Consequently  (9.6)  yields 
error  bounds  without  the  need  for  computing  with  vectors  in  R" . 

Theorem  A.  If  F  is  symmetric  with  respect  to  a  symmetric  positive  definite  A  then,  for  any 
y  €  R"\0  and  0e  R,  there  is  an  eigenvalue  X  of  F  satisfying 

lX-01  SIIFy-y0HA/llyllA  ■ 


Proof : 

llyllA  =  ll(F-0)->(F-0)yllAS 

ll(F  -  0)-‘|Ia  IIF  y  -  y  011^  =  IIF  y  -  y  011^  /  minv  I X  -  0 1  □ 

This  result  depends  strongly  on  the  existence  of  a  full  set  of  A-otthogonal  eigenvectors  for  F  and  is 
not  valid,  or  necessary  meaiungfiil,  when  A  is  indefinite.  Thus  we  cannot  use  (9.6)  just  as  it  is. 

For  indefinite  A  we  may  invoke  the  results  presented  in  [Kahan,  Paiiett,  and  Jiang,  1982]  and  adapt 
them  to  our  problem.  Two  residual  nonns  are  required.  Let 
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u:=Fy-y0  .  ^  F-ei?  ,  (9.7) 

where  y‘  x^O,  liyll  s  llxll » 1.  Now  we  are  obliged  to  woik  in  C”  rather  than  R” . 

Theorem  B.  \^th  the  notation  given  above  (6,  y.  x*)  is  an  eigentriple  of  F  -  E  for  some  E  satisfy¬ 
ing 

IIEII:Smax{  ilull.lfv'll )  . 

Note  that  Itv'  II  s  iivll.  In  other  words,  if  Hull  and  llvll  are  both  tiny  then  (6,  y,  )  are  an  eigentriple  for  an 
operator  almost  indistinguishable  from  the  given  one. 

In  our  case,  if  F  =  B"’  A  then  F*  =  A  B“*  and  we  have  from  (9.4), 

B“‘  A  y  -  y  0  =  (^y+,  /  (fy>,)  sfy)  • 

The  second  residual  is  obtained  easily  by  premultiplying  by  A. 

(A  B-’)  (A  y)  -  (A  y)  0  =  A  qj^i  (p,>,  /  co^.,,)  sfy )  . 

Let 


PO  )  =  (Py  +l/<0;+l)  IsO  )'  • 

Then  (0,  y,  A  y)  is  an  eigentriple  of  F  -  E  for  some  E  satisfying 


IIEII  ^  PO ) 


1  HA  q;>ill 
llyir  llAyll 


(9.8) 


(9.9) 


When  the  computation  has  proceeded  enough  that  I PO  )  I  very  small  then  first  order  perturbation 
theory  may  be  invoked  to  obtain  an  accurate  error  estimate  for  0  by  regarding  F  as  a  change  to  F  -  E.  If 
X  is  the  eigenvalue  of  F  closest  to  0  we  have,  in  general, 

ix-ei  ;=  -K7(IIEII^)  ,  (9.10) 

l^xl  Hyll  - llxll 

The  first  term  on  the  left  is  independent  of  E  and  is  the  condition  number  of  0  as  an  eigenvalue  of  F  -  E: 

cond(0)  =  llyir  llxll /  ly*  xl  .  (9.11) 


We  now  adapt  these  results  to  our  problem.  Note  that 

llyll^s  Is*  Qj  Qy  si  \t  s  \  =  j 
and 

lly'  A  yll  =  I  s'  Qy  A  Qy  si  S  I?  Qy  si 
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Now  we  can  obtain  an  error  bound.  As  IIEII  0.  using  (9.9) 


I  A,  -  0 1  ^  cond(e)  •  IIEII  +  O  (IIEIl^)  . 


Hyll  •  IIA  yll 
ly*  Ayl 


1  MA_q^l 

llyll  ■  IIA  yll  1 


+  <9(IIEI|2)  , 


—  max(  IIA  y  1 1  .IIA  q- jll  •  llyll }  +  O  (IIEII^)  . 
I?  Q;  si 


Is*  Qj  si 


IIAII  +  0(11EI|2) 


(9.12) 


The  attraction  of  the  last  inequality  is  that  the  dominant  term  can  be  calculated  at  step  j ,  without  recourse 
to  n -vectors,  provided  that  HAM  is  provided  along  with  the  subroutine  that  multiplies  vectors  by  A.  It  is 
useful  to  compare  (9.12)  with  (9.6)  and  Theorem  A. 


We  use^')  /  1^  Qj  sl  asa  provisional  error  estimate.  When  the  required  number  of  Ritz  values  6 
have  passed  this  test  then,  and  not  before,  the  Ritz  vector  y  may  be  computed.  At  that  point  the  more  pre¬ 
cise  factor 


max(  IIA  yll.  IIA  q^^jll  *  llyll ) 
may  be  computed  at  the  cost  of  forming  A  y. 

At  the  end  of  a  Lanczos  run  a  multiple  of  A  q^^.]  is  available  and  its  Euclidean  norm  can  be  com¬ 
puted  at  the  cost  one  dot  product  No  extra  call  on  A  is  necessary  for  that  term. 


10.  NUMERICAL  EXAMPLES 

In  this  section,  we  use  several  Test  Problems  to  assess  the  performance  of  the  proposed  algorithm  to 
extract  the  eigenpairs  of  damped  dynamic  systems.  The  mass,  damping,  and  stiffness  matrices  of  discre¬ 
tized  systems  are  computed  using  the  FEAP,  a  Finite  Element  Analysis  Program,  [Zienkiewicz  1977, 
Chapter  24].  The  results  reported  herein  are  obtained  using  a  VAX  Station  II/GPX  computer  system 
using  the  Ultrix  1.2  operating  system  and  the  f77  Fortran  compiler.  To  show  the  savings  resulting  from 
the  use  of  selective  orthogonalization  (SRO),  we  have  added  a  full  re-orthogonalization  (FRO)  option  to 
our  program  and  present  the  comparisons  in  Table  1. 

Test  problem  1  :  The  structure  is  modeled  as  a  cantilever  beam  with  a  lumped  translational 
viscous-damper  attached  at  the  tip.  The  beam  is  modeled  using  the  elementary  beam  theory  where  the 
geometrical  configuration  and  physical  properties  are  shown  in  Figure  2.  The  consistent  mass  is  used  to 
define  M.  The  damping  matrix  C  has  only  one  nonzero  element  representing  the  magnitude  c  of  the 
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lumped  damper.  The  cantilever  beam  is  divided  into  20  equal  elements  and  has  40  degrees  of  freedom. 
The  order  of  the  associated  (A,  B)  is  80. 

We  use  the  Lanczos  algorithm  with  the  FRO  scheme  to  solve  this  problem.  Figure  3  summarizes 
the  results  of  8  experiments.  Here,  we  call  a  Ritz  pair  good  if  the  pseudo  length  of  its  associated  residual 
vector  is  less  than  KT*.  This  criterion  ensures  that  a  good  Ritz  pair  approximates  the  wanted  eigenpair  to 
high  accuracy.  From  the  results  in  Figure  3,  we  see  that  the  first  few  eigenpairs  can  be  extracted  at  a 
fairiy  low  cost  compared  to  the  other  eigenpairs.  This  is  because  the  re-orthogonalization  cost  is  greater 
at  later  steps  in  the  Lanczos  algorithm. 

In  this  problem,  we  have  run  the  algorithm  to  compute  all  eigenvalues  in  order  to  test  the  robusmess 
of  the  computer  program  developed.  However  we  emj^iasize  that  the  algorithm  is  intended  only  for  par¬ 
tial  solution  of  a  large  eigenproblem.  After  the  80  steps,  we  see  that  the  pseudo  length  of  the  81th  Lanc¬ 
zos  vector  is  0.9x10"*^,  which  is  at  the  round-off  level,  implying  that  the  computed  Lanczos  vectors  have 
spanned  the  whole  solution  space  as  they  should  in  exact  arithmetic.  This  desirable  result  will  ensure  that 
all  the  Ritz  pairs  obtained  finom  the  solution  of  the  reduced  tri-diagonal  system  are  good  and  hence  are 
accurate  eigenpairs  of  the  system. 

Test  Problem  2  :  The  system  consists  of  two  beams  connected  by  a  hinge  with  a  rotational 
viscous-damper.  The  geometrical  configuration  and  physical  propenies  of  the  system  w  shown  in  Fig¬ 
ure  4.  The  consistent  mass  matrix  is  used  for  M.  The  damping  matrix  C  has  only  four  nonzero  elements, 
which  are  due  to  the  lumped  rotational  damper.  The  system  is  divided  into  40  equal  elements  and  has  83 
degrees  of  freedom.  The  associated  (A,  B)  is  of  order  166.  The  system  is  unrestrained  and  has  rigid  body 
modes,  so  we  apply  shift  to  the  (A,  B)  to  compute  the  eigenpairs  of  this  unrestrained  system.  The  Lanc¬ 
zos  algorithm  with  FRO  scheme  is  used  to  solve  this  problem.  Figure  5  summarizes  the  results  of  9 
experiments.  Similar  conclusions  as  in  the  first  test  problem  can  be  inferred  from  Figure  S.  The  pseudo 
length  of  the  167th  Lanczos  vector  is  0.1x10“*^,  which  again  exhibits  the  robustness  of  the  computer  pro¬ 
gram  developed. 

In  general,  the  starting  vector  for  the  Lanczos  algorithm  may  be  chosen  arbitrarily.  However,  if  the 
starting  vector  is  orthogonal  to  any  of  the  eigenvectors  of  (A,  B),  all  the  Lanczos  vectors  will  also  be 
orthogonal  to  these  eigenvectors.  In  practice,  round-off  errors  eventually  will  introduce  components 
along  these  eigenvectors;  however,  round-off  enters  slowly  and  the  convergence  to  these  eigenvectors  is 
deferred.  Therefore,  we  need  to  avoid  the  possibility  of  the  starting  vector  being  orthogonal  to  the  wanted 
eigenvectors  of  the  system.  Since  the  structural  system  in  this  test  problem  is  symmetric,  there  are  anti- 
sytmnetric  modes  as  well  as  symmetric  modes.  If  a  symmetric  starting  veaor  is  used,  such  as 
( 1, 1, . . . ,  1 ),  all  the  Lanczos  vectors  will  be  symmetric.  Accordingly,  all  the  anti-symmetric  modes  of 
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the  stnictuie  will  be  suppressed  by  this  biased  starting  vector.  To  obtain  all  the  required  lower  modes,  we 
must  avoid  choosing  either  a  symmetric  or  an  anti-symmetric  starting  vector.  This  is  a  strong  reason  for 
using  a  random  vector  to  start. 

Test  Problem  3  :  This  problem  is  a  three  dimensional  space  truss  system.  There  are  44  nodes  and 
the  4  erxl  tKxles  are  fully  restrained,  as  shown  in  Figure  6.  Thus,  there  are  120  degrees  of  freedom  and  the 
associated  (A,  B)  is  of  order  240.  All  truss  bars  have  the  same  density  and  Young’s  modulus  but  different 
damping,  as  shown  in  Figure  6,  resulting  in  a  nonproportionally  damped  system.  We  use  the  Lanczos 
algorithm  with  the  FRO  scheme  to  generate  60  Lanczos  vectors.  We  also  use  the  Lanczos  algorithm  with 
the  proposed  SRO  scheme  to  generate  60  Lanczos  vectors.  The  results  from  the  two  schemes  are  com¬ 
pared  in  Table  1.  The  SRO  scheme  is  shown  to  be  adequate  to  compute  the  desired  solution. 

Test  Problem  4  :  This  problem  is  a  larger  three  dimensional  space  truss  system.  There  are  300 
nodes  and  the  4  end  nodes  are  fully  restrained,  as  shown  in  Figure  7.  A  typical  cell  is  the  same  as  the  typ¬ 
ical  cell  in  the  test  problem  3.  There  are  888  degrees  of  freedom  and  the  order  of  the  associated  (A,  B)  is 
1776.  We  use  the  Lanczos  algorithm  with  the  FRO  scheme  to  generate  80  Lanczos  vectors.  We  also  use 
the  Lanczos  algorithm  with  the  SRO  scheme  to  generate  80  Lanczos  vectors.  The  results  from  the  two 
schemes  are  also  compared  in  Table  1. 

From  Table  1,  we  see  that  the  60  Ritz  pairs  obtained  provide  28  good  eigeipairs  for  test  problem  3 
and  the  80  Ritz  pairs  obtained  provide  40  good  eigenpairs  for  test  problem  4  for  both  FRO  and  SRO 
cases.  That  is,  approximately  two  Lanczos  vectors,  on  the  average,  are  required  to  capture  a  new  eigen¬ 
vector  for  these  two  large  problems.  This  implies  that  the  Krylov  subspace  generated  by  B~‘ A  and  a  ran¬ 
dom  vector  is  very  effective  in  approximating  the  least  dominant  eigenveaois  of  the  damped  dynamic 
systems  considered. 

By  maintaining  semi-orthogonality  between  the  Lanczos  vectors  with  the  SRO  scheme,  the  result¬ 
ing  Riu  values  are  as  accurate  as  those  obtained  with  the  FRO  scheme,  as  shown  in  Table  1.  But  a  great 
pan  of  the  re-orthogorudization  steps  can  be  eliminated  by  using  the  SRO  scheme  instead  of  the  FRO 
scheme.  That  is,  we  can  eliminate  partial  re-orthogonalization  effon  without  sacrificing  accuracy  of  the 
final  solution  when  solving  X  A  z  =  B  z  with  the  SRO  scheme.  This  is  in  agreement  with  the  case  of  solv¬ 
ing  o)^  M  w  s  K  w  by  a  standard  Lanczos  algorithm  with  the  SRO  scheme. 

To  assess  the  efficiency  of  the  Lanczos  algorithm,  the  lower  mode  solutions  of  the  above  four  test 
problems  are  also  computed  using  a  subspace  iteration  algorithm.  The  subspace  iteration  algorithm 
reported  in  [Chen  and  Taylor,  1986]  is  used  for  this  purpose.  The  subspace  dimension  is  determined  by 
min(2n,  n-t-8),  where  n  is  the  number  of  wanted  eigenpairs.  Table  2  compares  the  cost  of  the  Lanczos 
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Table  1  Results  of  Test  Problems  3  and  4 


item  problem 

Test  Problem  3 

Test  Problem  4 

FRO 

SRO 

FRO 

SRO 

number  of  Lanczos 

vectors  generated 

60 

60 

80 

80 

number  of 

re-orthogonalizations 

1770 

602 

3159 

1246 

CPU  time  spent  on 

generating  Lanczos  vectors 

40.3 

32.6 

536.6 

473.5 

CPU  time  spent  on  solving 

reduced  eigenproblem 

50.8 

51.7 

100.9 

101.2 

total  CPU  time  spent  on 

solving  the  whole  problem 

113.9 

106.7 

893.2 

830.9 

number  of  good 

Ritz  pairs  obtained 

28 

28 

40 

40 

algorithm  with  the  cost  of  the  subspace  iteration  algorithm.  It  is  apparent  that  the  Lanczos  algorithm  is 
considerably  more  efficient  than  the  subspace  iteration  algorithm  for  the  examples  considered.  However 
more  sophisticated  versions  of  subspace  iteration  might  perfonn  somewhat  better  than  ours  but  not 
enough  to  alter  the  striking  contrasts  in  Table  2. 

Breakdown  occurs  when  Uy  =  0  for  some  j .  The  algorithm  provides  a  bad  basis  if  there  are  any  co^ 
as  small  as  =  KT*.  In  Figure  8  we  plot  the  sign(fi}j)*  log  ^  against  j .  The  result  is  typical  for  our 

examples.  Quite  quickly  I  (Oj  I  drops  to  10"^  but  seems  to  stay  at  that  level  without  deteriorating  for  120 
steps.  We  have  no  explanation  of  this  phenomenon. 
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Table  2  Results  from  different  algorithms 


Test 

Lanezos  algorithm 

subspace  iteration  algorithm 

Problem 

No.  of  good 

Ritz  pairs 

CPU  time 

(second) 

No.  of  good 

Ritz  pairs 

CPU  time 

(second) 

1 

8 

8.5 

8 

40.2 

2 

20 

41.7 

16 

214.6 

3 

28 

113.9 

24 

1012.9 

4 

40 

893.2 

40 

20992.8 
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Figure  1  Outline  ofLanczos  inner  loop  with  operator  H"*A 


Step  1  : 

Operation  Count 

pick  a  starting  vector  r  (default  is  random) 

p  =  A  r 

2p(M)  +  p(C) 

solve  H  q  s  p 

v(K) 

P,=  lqT,l‘'2 

2n 

q<-q/Pi 

2n 

p  =  Aq 

2p(M)  +  p(C) 

©1  =  (q”^  p) 

solve  H  r  =  p 

v(K) 

«!  =  (r'*'  p) 

2n 

r  <-  r  -  q  (tti  /  CDi) 

2n 

oldp  =  A  r 

2p(M)  +  p(C) 

P2=  Ir'^rl^^ 

2n 

0)2  =  (r^  oldp)/ 1  r  1 

store  q  as  qi 

Loop:  Forj  =  2,3... 

Operation  Count 

oldq  4-  q 

oldp  < — >  p 

q  =  r  /  Py 

2n 

p  ^  p  /  Py 

2n 

solve  H  r  =  p 

v(K) 

(r^oldp  should  equal  Py) 

r «-  r-oldq  (py  /o)y_i) 

2n 

a,  <-  r'*'  p 

2n 

r<-r-q(ay  /  Wy) 

2n 

oldp  =  A  r 

2p(M)  +  p(C) 

py*,=  lr’'rl’'2 

2n 

0)y+j  =  (r^  oldp)/ 1  r  1 

(repeat  to  maintain  local  orthogonality) 

store  q  as  qy 

Figure  5  Results  of  Test  Problem  2 


Figure  6  Test  Problem  3 


Figure  7  Test  Problem  4 


sign  (a>i)*  log 


